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STRACT 

his paper we propose several strategies for the 
exact computation of the determinant of a ratio- 
n«F~matrix. First, we use the Chinese Remain- 
dering Theorem and the rational reconstruction to 
re&ever the rational determinant from its modular 
images. Then we show a preconditioning for the 
determinant which allows us to skip the rational 
reconstruction process and reconstruct an integer 
resul t. We compare those approaches with matrix 
preconditioning which allow us to treat integer in- 
stead of rational matrices. This allows us to intro- 
ducc integer determinant algorithms to the ratio- 
nafcdeterminant problem. In particular, we discuss 
tr^ applicability of the adaptive determinant algo- 
rithm of [9] and compare it with the integer Chinese 
Remaindering scheme. We present an analysis of 
the, complexity of the strategies and evaluate their 
Ojrpc rimental performance on numerous examples, 
fois experience allows us to develop an adaptive 
fftT^tegy which would choose the best solution at 
frBT^run time, depending on matrix properties. All 
strategies have been implemented in LinBox linear 
algebra library. 

INTRODUCTION 

TJje determinant computation is one of the core 
problems in linear algebra. To our knowledge, the 
problem of the exact computation of the determi- 
nanjt of a rational matrix (i.e a matrix with ratio- 
nal entries) has not yet been widely studied. In 



general, exact algorithms can be used everywhere 
where large precision is required. For example, the 
determinant can be too close to or ±00 and thus 
cannot be computed by floating point precision al- 
gorithms. In the case of ill-conditioned matrices 
symbolic methods can be preferred as rounding er- 
rors can spoil the computation. It can also be in- 
teresting to compare the use of decimal and contin- 
ued fractions approximations of the entries of real- 
valued matrices. Continued fractions are the best 
approximants with small denominators, see [12, Ch. 
4]. In this paper, we will try to face the question 
of how efficient an exact determinant computation 
can be in both cases. 

LinBox library [7] implements exact algorithms for 
the determinant computation in the case of mod- 
ular and integer domains. By using fast modular 
routines [6, 8] it can offer solutions an order of mag- 
nitude faster than other existing implementations 
[9] . We apply these procedures to the computation 
of the determinant of a rational matrix. 

Rational field arithmetics is implemented in CMP 
[21] and Givaro [22] libraries. In general, rational 
numbers are difficult to treat from the exact com- 
putation point of view. Mainly, the size of the nu- 
merator and denominator can increase very quickly 
with every addition and multiplication. When we 
add or multiply two fractions with numerators and 
denominators bounded by M, the numerator and 
denominator of the result are bounded by 0(M 2 ). 
Moreover, one addition requires 3, and one mul- 
tiplication requires 2 integer products, as well as 
a gcd computation. Therefore, the cost of an ex- 
act matrix-vector or matrix-matrix product can be 
prohibitive in practice. This prohibits the use of 
the rational field Q in most exact linear algebra 
algorithms which rely on matrix-matrix or matrix- 
vector products. 



However, the cost of computing a modular image 
of a rational number £, where a, b are of moderate 
size, should be comparable with the cost of com- 
puting a modular image of a large integer number. 
This allows us to compute a modular image of a ra- 
tional matrix at a reasonable cost and thus enables 
us to use modular procedures. 

To compute the determinant of a rational matrix 
A = [j^-],bij > the problem of matrix storage 
has to be considered. First, we can store the entries 
of A as rational numbers. Furthermore, one could 
store the common denominator D(A) of all entries 
of A and an integer matrix A' given by the formula 
A = gm A'. This approach can be useful in the 
case when the entries of A are decimal fractions and 
D(A) can be set to a power of 10. But if we only 
assume that the values \cnj\,bij are less than M, 
both D(A) and p'|| are bounded by 0(M™ 2 ). Still, 
we may store the common denominator for each row 
(column) separately. Then the integer vectors Ai 
are given by the equation Ai = -k-Ai, where Ai 
is the matrix row (column) and Di is the common 
denominator of its entries. Vectors Ai form matrix 
A, the norm of which is bounded by 0(M"). The 
product irDi gives a more accurate approximation 
of the denominator D(det(A)) than D(A) n . 

The purpose of this paper is to propose the strate- 
gies to compute the denominator of a rational ma- 
trix. All approaches are based on modular com- 
putation. Depending on the matrix storage deter- 
minant and/or matrix preconditioning is proposed. 
The resulting algorithms can use the rational re- 
construction [12, Ch.5] and/or existing integer de- 
terminant algorithms. 

The rest of the paper is organized as follows. In sec- 
tion 2 we give a short description of the existing al- 
gorithms for the rational reconstruction and the in- 
teger determinant problem. In section 3 we present 
the main result i.e. two preconditioning strategies 
and four new algorithms to compute the rational 
determinant. The cost of the algorithms can be 
described in terms of the number of modular im- 
ages of A and modular determinant computations 
needed. Depending on the strategy, the cost of the 
rational reconstruction or p-adic lifting is taken into 
account. In section 4 we discuss the cost of com- 
puting a modular image of a matrix and the overall 
cost of the algorithms. In section 5 we present the 
experimental results and discuss the best choice of 
the strategy in practice. We conclude the paper by 
proposing some mixed solutions in section 6. 



2. EXISTING ALGORITHMS 

The aim of this section is to introduce the algo- 
rithms that will be used later in section 3. In sub- 
section 2.1 we give a short description of the ratio- 
nal reconstruction procedure. On the example of p- 
adic system solving [2], we present the application 
of this procedure to the computation of a rational 
solution. We show how to change the procedure in 
the case of early terminated reconstruction [18] and 
give the complexity estimation in this case. Then 
in subsection 2.2 we present the classical CRA al- 
gorithm for the determinant and its modifications 
by [1] and [9]. 

2.1 Rational reconstruction and its ap- 
plication 

A modular image of a rational number ^ mod M 
can be computed by taking the modular images of 
a and b and applying the modular division. This 
fact can be written as 

— — u mod M <4> a — bu mod M. 
b 

It should be noticed is that the opposite procedure 
can also be performed. One can reconstruct the 
fraction | where gcd(a, 6) = 1,6 > from it mod- 
ular image u. The solution is usually not unique 
but when we additionally require that |a| < -y, 
b < 77, then there exists at most one solution, see 
[12, Ch.5]. 

The solution to the rational reconstruction prob- 
lem can be computed by applying the extended 
Euclidean algorithm EEA which searches for the 
gcd of M and u. The procedure Ratrec(a,b,u,M, 
N, D) takes as the input modulus M, u € Z and 
the bounds N and D, and returns a fraction f = 
umodM such that |a| < N,b < D or FAIL if no 
such solution exists. The worst case complexity 
of Ratrec is thus the same as for the EEA algo- 
rithm i.e. 0(log 2 (M)) for the classical algorithm 
and C>(log(M)log(log(M))) for the fast Euclidean 
algorithm, see [12, Ch.ll]. We will use the nota- 
tion EEA(M) for the complexity of the Extended 
Euclidean Algorithm with entries bounded by M. 

In many application, the cost of rational recon- 
struction is usually small compared with the cost of 
computing u and M. The general scheme is to re- 
cursively compute Wfe, M fc , where M k = p\pi ■ ■ ■ Ph 
or M fe = p k until M k > 2ND and then to apply the 
rational reconstruction. The complexity of the pro- 
cedure depends on the number k of steps, which can 
be quite large Reducing the number of steps can be 
the easiest way to enhance the performance of the 



algorithm. 



9] for a detailed description. 



This can be seen on the example of the Dixon al- 
gorithm [2] to solve a linear system Ax = b of in- 
teger equations. Let N, D be the bound for the 
numerator and denominator of a;. In the classi- 
cal approach we compute the p-adic approxima- 
tion in k > log(N) + log(D) + 1 steps and then 
reconstruct the result, which gives the complex- 
ity O (m 3 (log(m) + log(P||)) 2 ) when we use the 
bound of Hadamard for D, N and assume b £ O(l). 
See [15] for a detailed complexity study. In fact, the 
number of entries in x which we need to reconstruct 
can often be reduced, see [7]. 

One should however notice, that the bounds N and 
D can be much bigger than the actual result. The 
idea is therefore to apply the rational reconstruc- 
tion periodically and check the solution for correct- 
ness. If Mfe = p k is the modulus in the current step, 

the method of Wang [18] prompts us to set 
as the current bound for numerators and the de- 
nominator in Ratrec. The algorithm is guaranteed 
to return the result if M k > 2ma,x(N(x) 2 , D(x) 2 ), 
where N(x),D(x) are the values of the numera- 
tor and the denominator. In the opposite case, 

Ratrec(a, b, u, Mk, \J -^S \J should fail with large 
probability. If we apply Wang's idea to the p-adic 
lifting we can reduce the number of steps to k = 
2 log p (max(iV(a:), D(x))) + 1 and the complexity be- 
comes O (m u + m 2 fclog(m|| J 4||)) Current work on 
this field focus on further reducing the number of 
steps in the case when N(x) <C D(x) or D(x) <g; 
N(x). A purely heuristic idea is to use the bounds 

\ft¥%> instead of For other ap- 

proaches, see [14, 17]. 

2.2 Integer Determinant Algorithms 

For an integer matrix A one has several alterna- 
tives to compute the determinant. The classical ap- 
proach is to use Chinese Remaindering Algorithm 
(CRA) to reconstruct the value from sufficiently 
many modular images. The modular determinant 
is computed by LU factorization in the time 0(n"), 
where n is the matrix dimension. Each step of the 
algorithm consist of computation mod pi and a re- 
construction of the determinant mod p\ ■ ■ ■ pi by the 
Chinese Remaindering Theorem. The computation 
is stopped when the early termination (ET) con- 
dition is fulfilled i.e. the reconstructed result rests 
the same for several iterations. The algorithm is 
Monte Carlo type, where the probability of success 
is controlled by the number of repetitions. See [4, 



A mixture of CRA loop and Dixon p-adic lifting is 
used to compute the integer determinant in [1] and 
in the hybrid algorithm of [9]. The principle is to 
reduce the value reconstructed by CRA algorithm 
by computing a large fraction of the determinant. 
By solving several linear systems we can compute 
some largest invariant factors s m , . . . s m -i- Their 
product 7r is potentially a large part of the determi- 
nant. An early terminated CRA loop which recon- 
structs det(A)/n mod poPi ■ ■ ■ Pi usually requires 
only a few modular determinant computations. In- 
formally, the algorithm can be described as follows. 

1. For i — to k do 

(a) Solve Axi = bi by Dixon p-adic lifting to 
find 

(b) TV — S m ■ ■ ■ S m — j, 

(c) Run CRA for several iterations to deter- 
mine det(A)/iv; 

(d) if ET break; 

2. Run another determinant algorithm to get the 
result; 

Here, k should not exceed the e xpected number 
of invariant factors which is 0(y / Iog(n)) see [9]. 
The expected complexity of the hybrid determi- 
nant algorithm [9] for random dense matrices is 
O (n 3 log 2 - 5 (n||A||)). In the worst case (step 2) we 
can choose between the CRA algorithm and the al- 
gorithms of [16, 10, 13]. In fact, in the expected 
case we do not need to run this step. The experi- 
ment proved that thanks the adaptive solutions this 
algorithm performs better than other implementa- 
tion for a larger group of matrices. 

3. RATIONAL DETERMINANT AL- 
GORITHMS 

The algorithms to compute the rational determi- 
nant are based on the ideas described in section 
2. We present four main strategies to compute the 
rational determinant. They all use CRA which al- 
lows us to compute the determinant of the matrix 
modulo a product pi - ■ - pk of primes. Then the 
first variant uses the rational reconstruction to ob- 
tain the rational result. In order to make use of 
Early Termination condition we have to precondi- 
tion the determinant to obtain its integer multipli- 
cation. Preconditioning of the matrix allows us to 
use the integers determinant algorithms. The ap- 
plication of two determinant algorithms is studied 



here. The common requirements for all algorithm 
are shown in 3.0. The algorithms are Monte Carlo 
type due to the early termination used. 



Requirements 3.0 

Require: A - an m x m rational matrix; 

Require: Di, i = 1 . . . m - the common denomina- 
tor of the entries of the ith row (column); 

Require: N,D - the bounds for the numerator and 
the denominator of det(A), D — ivDi; 

Require: A set P of random primes; 

Require: Ratrec(a, b, u, M, N,D) - a procedure 
which reconstructs | = u mod M,a < N,b < D 
or returns FAIL. 

Ensure: det(A) - the determinant of the matrix. 



Algorithm 3.1 RatLU 
1: i = 0,fc = 0,n = 0, d = 1,M = l,u = 0; 
2: repeat 

3: Get pi from P; 

4: Compute Ai — A mod pv, 

5: Compute Ui = det(Ai); 

6: Reconstruct u = det(A)modMp; using 

M,u,Ui,pi, M = Mp r ; 
7: if i = fc 2 then 

8: s = Ratrec(n, d, u, M, y^);++fc; 

9: if s 7^ FAIL then return n, d; end if 
10: end if 
11: until M > 2ND 
12: status = Ratrec(n, d, u, M, N, D); 
13: if status / FAIL then return n, d; end if 



The effectiveness of our methods depends heavily 
on the number of modular determinants computed 
and thus on the bound N and D for the numerator 
and the denominator of the determinant. One can 
compute D as the product of 1cm of all denomina- 
tors in a row (or a column). Then N can be com- 
puted as D ■ H, where H is the Hadamard bound 
for matrix A. One should notice that the bounds 
can be largely overestimated. Thus, we proposed 
output-dependant approach which allows us to re- 
duce the number of iteration. 

The first idea is to employ the CRA scheme and 
compute the determinant for the modular images 
of a rational matrix. In the case when the determi- 
nant is rational, early termination condition never 
holds. Instead, we have to compute the bounds 
D and N for the denominator and numerator of 
the determinant. As soon as the product of primes 
M = pi ■ ■ ■ pk overcomes 2ND we can apply ratio- 
nal reconstruction and reconstruct the determinant 
from the modular image. We can also use an output 
dependent rational reconstruction as described in 
section 2.1. This strategy is presented as algorithm 
RatLU. An early termination in the rational case 
would required applying the rational reconstruction 

from time to time with the bounds N = D — \J ^Jl 
and wait for the result to re-occur. This leads to so- 
lution when M > 2 max{n 2 , d 2 }, where n, d are the 
numerator and denominator of the determinant. 

The second method can use the denominator bound D 
to make the CRA loop look for an integer value. 
Again, we compute the modular image of a ratio- 
nal matrix A but this time we call CRA to look for 
D x &et(A) which is integer. Now the classic ET 
condition can be used and the result is obtained as 
soon as M > n-j. The effectiveness of this method 



depends therefore on the exactness of denominator 
bound D. Experimental results show that it is suf- 
ficient in practice, see sec. 5 table 2. This strategy 
is presented as algorithm PrecDetLU. 



Algorithm 3.2 PrecDetLU 
1: i = 0;M = l;u = 0; 
2: repeat 

3: ++i;Get pi from P; 

4: Compute Ai — A mod p» ; 

5: Compute m — D ■ det(Ai); 

6: reconstruct u — D ■ det(A)modMp; using 

M,u,u l ,p i , M = M-pi 
7: if ET holds then return " m , .? m ; 

gcd{u,D) ' gcd(u,D) ' 

end if 

8: until M > 2ND 

9: return gcd(UjD) , gcd ( u ,D) ' 



The last two strategies require an integer matrix A 
which can be obtained by preconditioning the ra- 
tional matrix A. In order to obtain an integer ma- 
trix, the easiest way would be to take matrix A' = 
D(A)A, where D(A) is the common denominator 
of all entries. In the general case, where the entries 
of A are fractions -r 11 with numerator and denom- 

bij 

inator bounded by ||A||, this is not the best choice 

2 

as the size of D(A) can be as large as 0(||A|| m ). 
This causes log(||A'||) to be 0(m 2 ). Moreover, the 
denominator approximation is D(A) m in this case, 
which is 0(m 3 ) in size. We have already defined 
a tighter bound for the denominator of det(A) by 
irDi, which is 0(m 2 ) in size. Now, if we want to 
use the integer matrix A then we can precondi- 
tion A by taking A = Adiag(Di), where Di are 
the common denominators of the rows (or A = 
diag(Z);)A, where Di are the common denomina- 



tors of the columns). For the preconditioned ma- 
trix A all integer determinant algorithms can be 
applied. In particular the hybrid determinant algo- 
rithm of [9] can be used. The drawback of this ap- 
proach is the size of the coefficients of A compared 
to A, see section 5 table 1. This forced us to use 
early terminated rational reconstruction for system 
solving in the Dixon p-adic lifting algorithm. The 
strategies that use the CRA loop or the hybrid algo- 
rithm are presented as algorithms PrecMatLU and 
PrecMatDixon respectively. 



Algorithm 3.3 PrecMatLU 
1: i = 0;M = l;u = 0; 

2: Compute A = Adiag(A) (or diag(A)A) 

3: repeat 

4: Get pi from P; 

5: Compute Ai = A mod p» ; 

6: Compute m = det(Aj); 

7: reconstruct u = det(A) mod Mpi using 

M,u,Ui,pi, M = M ■ pi 
8: if ET holds then return " m , .. g m ; 

end if 
9: until M > 2ND 
10: return gcd( " D) , gcd(lt C) ; 





Alj 


gorithm 3.4 PrecMatDixon 


1: 

2: 
3: 


Compute A = Adiag(A) (or diag(A)A); 
Compute u = det(A) by HybridDet [9]; 
return u D 



4. COMPLEXITY ANALYSIS 

In this section we study the complexity of the al- 
gorithms presented in section 3. In subsection 4.1 
we present the analysis of the general case, where 
we assume that the entries of the matrix are frac- 
tions with numerators and denominators bounded 
by ||A||. Then, in subsection 4.2, we will focus on 
two special cases i.e. matrices of decimal fractions 
and Hilbert matrices. 

The complexity of the strategies described in sec- 
tion 3 depends on the number of iterations required 
by the while loop of CRA. Then, depending on the 
strategy, we have to include the cost of computing 
the homomorphic image of the matrix, the cost of 
the rational reconstruction or the cost of p-adic lift- 
ing. If we use the early termination condition, the 
number of steps required for the computation of 
det(A) depends on the values: m - the size of the 
matrix, n, d - the real values of the numerator and 
denominator of det(A) and D - the bound for the 
denominator. The cost of homomorphic imaging 



depends on the maximum norm of the matrix i.e. 
\\A\\ = max{ || a,ij ||, bij] and ||A||. 

4.1 General case 

We start this section by the analysis of the ratio- 
nal homomorphic imaging schemes. We have the 
following lemma. 

Lemma 4.1. Let p be a word-size prime. Then 
the complexity of computing the modular image at p 
for a rational matrix A is 0(m 2 (log(|| A||))+EEA(p)) 
word operations. 

Proof. For a matrix without a pattern we com- 
pute an image for all m entries. For a rational 
fraction the cost is 0(log(||A||)) for the computa- 
tion of the modular image of the numerator and de- 
nominator and EEA(p) = 0(log(p) log(log(p))) for 
the modular inverse computation by fast extended 
Euclidean algorithm. Therefore for a word-size ||A|| 
the cost of computing the image is O(l) yet impor- 
tant, due to the constant for computing the inverse 
of an element mod p. □ 

For the integer case, the cost is log(||A||)). We 
can notice that log(||A||) can be 0(m log(||A||)) in 
the worst case, so the complexity of homomorphic 
imaging in terms of m is 0(m 2 ) for the rational 
and 0(m 3 ) in the integer case. But if ||A|| < p the 
cost of imaging for one element is 1. Thus, if both 
|| A|| and || A|| are less than p, the complexity of the 
homomorphic imaging becomes m EEA(p) for the 
rational and m 2 for the integer case. In this case, 
it is better to use integer imaging. On the other 
hand, if matrix A is structured, for example it is 
Hankel-type, we have the complexity mEEA(p) for 
rational imaging. Due to the preconditioning, we 
loose the structure pattern for A and the complex- 
ity of integer imaging rests without change. Finally 
we notice, that for sparse matrices with Q elements, 
we can take O instead of m 2 in the complexity for- 
mula. 

Putting it together we have the following theorem. 

Theorem 4.2. The worst case complexities of the 
strategies for computing the determinant of a ratio- 
nal matrix A of size m are 

1. O (k(m 2 log(|| A||) +m"))+0~ (kV~k) for Rati U, 
where 0~ hides some log(fc) factors; 



2. O (log(f n)(m 2 log(P||) +m w )) for PrecDetLU; 

3. O (log(f n)(m 2 log(||i||) + m w )j for PrecMatLU; 

4. 0~( a; (m 2 (log(m)+log(||i||))+m i -*)+0(log(f £ + 
l)(m 2 log(||A||)+m w )) for PrecMatDixon, where 
s-m = s m (A) and x € m(log(m||j4||||6||) is the 
size of solution to Ax — b. 

Here A is equal to Adi&g(Di) as in section 3; n, d 
are the numerator and denominator of det(A) and 
k — 0(max(log(n), log(d))). 

PROOF. The complexities can be obtained by a 
careful examination of the number of CRA steps. 
The result for alg. RatLU takes into account the 
cost of the rational reconstruction which is per- 
formed at most 0(Vk) times. In alg. PrecMat- 
Dixon we introduce x to estimate the cost of early 
terminated p-adic lifting. The size of x can gen- 
erally vary depending on the choice of b but is 
0{m log(m|| A\\ \\b\\)) in the worst case. To further 
evaluate the worst case complexity of alg. PrecMat- 
Dixon we assumed that HybndDet continues to use 
CRA loop in the worst case. Thus the number of 
iterations 0(log(-2-jp-)) and the complexity. □ 

Special care should be taken if we consider the use 
of alg. PrecMatDixon. As \\A\\ can potentially be 
0~(m) in size and with a pessimistic bound on 
x, its worst case complexity can be 0~(log(m 4 )), 
which is worse than for the CRA computation. Nev- 
ertheless, the gain of computing s m can be impor- 
tant, as it is the case in the HybridDet algorithm, 
see [9]. 

4.2 Complexity in the special cases 

By the precedent remarks it should be visible, that 
the analysis of the strategies should be divided into 
two main cases. One would consist of the matri- 
ces, whose entries are given by decimal fraction, 
or more generally, where the common denomina- 
tor of all entries, the common denominator of the 
rows and the norm of A are of the same order i.e. 
D(A) = 0(A) = O(PII)- In the other case ma- 
trix entries are given as fractions with different de- 
nominators. We will study the complexity of the 
algorithms on the example of Hilbert matrices. 

In the case of matrices of decimal fractions let us 
further assume that || A\\ is O(l). This would be the 
case of numerous ill-conditioned matrices emerging 



from different applications in science and engineer- 
ing. In order to better describe the differences be- 
tween the algorithms, we include the cost of EE A 
when it is relevant. The theorem is a straightfor- 
ward consequence of theorem 4.2. 

Theorem 4.3. The complexities of the strategies 
in the case when \\A\\ = 0(||A||) = O(l) are: 

1. C>~ (k{m 2 EEA(p) + m w + kVkj for alg. RatLU; 

2. 0~ (log(-jn)(m 2 EEA(p) + m")) for alg. PrecDetLU; 

3. 0~ (log(f n)(m 2 + m u )) for alg. PrecMatLU; 

4. 0~ (x(m 2 log(m) + m:r5)) + log(f ^)(m 2 +m")) 
for alg. PrecMatDixon. 

where k,x are as in theorem 4-2. 

The analysis suggests that the algorithm PrecMatLU should 
be better than PrecDetLU (see 4.1 for the homo- 
morphic image complexity) . The performance anal- 
ysis in section 5 confirms this observation. Further- 
more, as long as the Smith form of A is simple, we 
encourage the use of strategy PrecMatDixon. In 
particular, we can establish an equivalence between 
matrices A of random decimal fractions with e dec- 
imal places taken randomly an uniformly from the 
interval [0,1] and matrices A, \\A\\ < 10 e . This 
allows us to use the expected complexity of the hy- 
brid algorithms of [9] as the expected complexity 
of the rational determinant computation by alg. 
PrecMatDixon. Also, the preconditioning should 
be used instead of strategy RatLU. For more de- 
tails see section 5. 

The other group consists of matrices with ratio- 
nal entries given by fractions with very different 
denominators. As a model case we can consider 
Hilbert matrices. Hilbert matrices are the matrices 
of the form H m = [ j+J 1 _ 1 ]i,j=i..m. They are bench- 
marks examples for many numerical methods. The 
formula for the determinant of a Hilbert matrix is 
well known and is given by the equation 

5 ^- J .n = ?(a + 1»( 2 *^ 

Theorem 4.4. The complexities for rational de- 
terminant strategies in the case of Hilbert matrices 
are 



1. O [m 2 \og{m){m w + m^/log(m))j for alg. Rath 

2. O {m w+2 log(m)) /or aZ 5 . PrecDetLU; 

3. O (m 5 )log(m)) /or alg. PrecMatLU; 

4- 0(s m m 3 log 2 (m) + m 5 log(m)) /or a/g. Prec- 
MatDixon. 

PROOF. One should notice that log( det ^ j ) is 
0(m 2 log(m)). The size of entries of H m is log(||_ff m ||) 
0(log(m)) andlog(||/f m ||) = 0(mlog(m)). □ 

In the case of Hilbert matrices algorithm PrecDetLU - 
has the best time complexity and also performed 
best in the experiments, see section 5. Since the 
numerator is equal to 1, we only have to recover 
the size of the over-approximation. Experimental 
results show, that its size is equal to about 8% of 
the denominator size. Therefore, alg. PrecDetLU, 
PrecMatLU perform about 25 times less iterations 
than RatLU. As for the algorithm PrecMatDixon, 
the study of the Smith form of H m has revealed that 
it is quite complex, with about 2y / m nontrivial fac- 
tors and the size \og(s m (H m )) equal 0(m). Thus, 
it is not worth computing PrecMatDixon due to 
the high cost of the algorithm and poor gain. 

5. PERFORMANCE COMPARISON 

In this section we present the experimental results 
for four strategies from section 3. We have tested 
the performance of four strategies on three matrix 
sets: random, ill-conditioned and Hilbert matrices. 

We generated the random matrices using Matlab 
procedure rand. The entries of the matrices are 
decimal fractions with 6 decimal places chosen ran- 
domly from the interval [0, 1]. The determinant of 
the resulting matrices is large in the absolute value. 
The result of the numerical procedure of Matlab is 
±co. 

Ill-conditioned matrices have been chosen from the 
Matrix Market [20] Harwell-Boeing collection. We 
chose three sets: Grenoble, Astroph and Bcsstruc3. 
Grenoble set represents the results of the simula- 
tion of computer systems. The sizes of the matri- 
ces varies from 115 to 1107 and the condition num- 
bers range from 1.5 • 10 2 in the case of the smallest 
matrix to 9.7 • 10 7 for the biggest. The decimal 
precision of the entries depends on the matrix and 
ranges from 1 to 5 decimal places. The determi- 
nants are close to 0. For these matrices, Matlab 
procedure det computes the result correctly up to 



the 5th decimal place. Since matrix entries seem to 
be represented as rounded expansions of rational 
numbers, we computed the determinant of the ma- 
trices " as is" and then we took continued fractions 
approximants of the entries with the same precision 
as the decimal fractions. 

Astroph set describes the process of nonlinear ra- 
diative transfer and statistical equilibrium in as- 
trophysics. The condition number is 3.6 • 10 17 for 
the small 180 x 180 matrix and 1.7 • 10 14 for the 
765 x 765 one. The result of Matlab computation 
is — co. Bcsstruc3 gives dynamic analyses in struc- 
tural engineering. All matrices are symmetric. The 
condition number is about 10 11 for matrices 19 and 
20 and 10 5 for matrix 22. The result of Matlab 
computation is co. 

We split the analysis of the performance of the al- 
gorithms in three phases. First, we will consider 
the cost of rational-modular vs. integer-modular 
imaging and compare it with the results for ||^4|| 
and \\A\\. Then we will take a look on the nu- 
merator and denominator approximations D and N 
computed by our algorithms. Finally, we give the 
timings for all strategies and compare their perfor- 
mance. 

As we can see in table 1, the time of computing an 
integer image can be several times shorter than for 
the rational image provided that the size of precon- 
ditioned matrix is still small. This is not the case 
for Hilbert matrices of dimension > 250 , when the 
time of rational image computation is better. Fur- 
thermore, for structured matrices, like Hilbert, we 
can reduce the number of images computed. For 
a Hankel-type matrix, there are only 2n — 1 im- 
ages to compute, which makes the cost of imaging 
negligible. 

The performance of the algorithms depends on the 
accuracy of denominator approximation used. For 
the bound D = nDi, the resulting size of the over- 
approximation is shown in table 2, column 4. In 
algorithm PrecMatDixon we additionally approx- 
imate the numerator by computing s m (A). In this 
case we are interested in the value App(N) = s m (A) 
and -j which we compute instead of the nu- 

merator. As we can see in the table, the quality 
of the approximation of the denominator depends 
on the matrix and ranges from 1-2% in the case 
of sparse matrices in the Grenoble set, to 80% for 
Bccstk matrices. For Hilbert matrices the approxi- 
mation is quite efficient, the over-approximation is 
always less than 10%. Table 3 shows that despite 
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Table 1: Comparison of the times (in sec- 
onds) for homomorphic imaging are given in 
columns Ratlm (for rational) and Intlm (for 
integer). The ratio of the timings is given in 
column 3. Last two columns give the size of 
entries for A and A. Matrix size is included 
in its name. 



the size of the over-approximation, preconditioning 
allow us to gain enough to beat the naive RatLU al- 
gorithm. If the size of ||A|| is small, as is the case 
for sparse matrices, we can compute s m (A) at a 
relatively low cost and efficiently approximate the 
numerator. 

The timings for all algorithms are shown in table 
3. The results for Hilbert matrices agree with the 
complexity estimation in Thm. 4.4. Note that alg. 
PrecMatDixon is usually the best for the matrices 
from MatrixMarket collection. 

For the Grenoble set, the approximation by contin- 
ued fractions allowed quite well, in our opinion, to 
reconstruct the orginal rational matrix connected 
to the problem. Despite the difference in proper- 
ties, the running times for the decimal and con- 
tinued fractions variants were simmilar. However, 
although the matrices were close in the maximum 
norm, the determinants ratio reached as much as 2 
in the case of grenoblell07. 

In figure 5 we present the results of the determinant 
computation for Hilbert matrices. We compare the 
timings for algorithm RatLU, PrecDetLU, Prec- 
MatLU, PrecMatDixon, and the Maple LinearAlge- 
bra::Determinant algorithm with method=rational. 
The best performance is observed for a variant of 
algorithm PrecDetLU which takes into account the 
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Table 2: The size of the numerator n and de- 
nominator d of &ct(A), the size of the denom- 
inator over-approximation D/d computed by 
PrecDetLU and PrecMatLU; the numera- 
tor approximation App(n) obtained as s m in 
PrecMatDixon, and the size of the part re- 
maining to compute. s m depends on n and 
the over-approximation D/d. 
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Table 3: Timing comparison for 4 rational 
determinant strategies. All times in seconds. 
Best times in bold. 
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Figure 1: Comparison of the timings for 
the exact computation of the rational de- 
terminant of Hilbert matrices. The results 
for algorithms RatLU, PrecDetLU, Prec- 
MatLU and PrecMatDixon implemented in 
LinBox and Maple Determinant procedure 
are shown. Algorithm PrecDetLU is used 
in the classic and symmetric variant, which 
takes into account the Hankel structure of 
the matrix. All times in seconds. 



1. Compute D = irDi, A; set TV = 1; 

2. If log p (||i|| <C) compute N = s m {\\A\\) - see 
alg. PrecMatDixon 

3. Compute the modular image of the rational 
matrix A and integer matrix A, determine 
whether to use PrecDetLU or PrecMatLU based 
on the timings. 

4. Run the ET CRA loop for § ■ det(A) using 
PrecDetLU or PrecMatLU. 

5. From time to time check by rational recon- 
struction the early termination condition on 
det(A) - see RatLU. 

This algorithm can be further developed to com- 
pute other invariant factors as in alg. PrecMat- 
Dixon if relevant. Notice, that the cost of intro- 
ducing solution RatLU to the adaptive algorithm 
is virtually that of rational reconstruction. 



Hankel structure of the matrix. 



6. CONCLUSIONS 

It this paper we have presented four strategies for 
exact computation of the determinant of a ratio- 
nal matrix. We have evaluated the performance of 
these algorithms on several sets of matrices. The 
performance of the algorithms suggests that there 
exists a clear division between the matrices given as 
a rational approximation (by decimal fractions) of 
real valued matrices and the matrices with a great 
diversity of the denominators of the entries. For 
the first case, matrix preconditioning which leads 
to a integer matrix is proposed, which allows us 
to use integer determinant algorithms, see solution 
PrecMatDixon. For the second case, determinant 
preconditioning is preferred, which does not lead 
to matrix coefficient blow-up. In general, precon- 
ditioning proved more useful than rational recon- 
struction tools, although better early termination 
methods where the modulus M is linear in the size 
of the output n and d can bring a change, see [14, 
17]. 

An adaptive solution should be able to choose the 
best storage method and homomorphic imaging scheme, 
and work independently of the determinant over- 
approximation. 

We propose the following solution, which incorpo- 
rates the elements of all algorithms 



Further work can include intertwining algorithms 
RatLU and PrecDetLU to include the use of less 
exact determinant preconditioners, which potentially 
are not a multiple of d. The aim would to re- 
duce a factor of the denominator by precondition- 
ing and reconstruct the remaining part by rational 
reconstruction. The strategy should be effective, 
if the over-approximation caused by precondition- 
ing is reduced but a large fraction of the denomi- 
nator is obtained at the same time. For example, 
D — nDi/ gcd(Di) could be considered. 

Further work can then focus on the implementation 
of the solution in the case of sparse matrices and 
on the parallelization of the algorithms. 

In this paper we have considered the case of dense 
matrices in the analysis of the complexity of the 
strategies as well as in the implementation. How- 
ever, sparse matrix counterparts of the algorithms 
can also be used. For the modular determinant 
computation one could used the algorithm of Wiede- 
mann [19] that computes the determinant by find- 
ing the characteristic polynomial of the matrix. In 
alg. PrecMatDixon the sparse solver of [11] can be 
used. 

The strategies described in this paper contain ele- 
ments that allow parallelization. This concerns in 
particular the CRA loop, where several iterations 
can be performed at the same time, see [3]. The 
question of an optimally distributed early termina- 
tion in the case of integer Chinese reconstruction 



(alg. PrecDetLU, PrecMatLU, PrecMatDixon) as 
well as the rational reconstruction (alg. RatLU) 
has not yet been addressed. For a parallel p-adic 
lifting for alg. PrecMatDixon, see [5]. 

In this paper we have developed and compared four 
strategies to compute the rational determinant of 
a matrix. We have proposed two preconditioning 
methods that allow us to transfer the problem from 
rational to integer domain. We believe that the ap- 
proach described in this article can also be applied 
in other problems of exact computation in ratio- 
nal numbers such as rank computation or system 
solving. 
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